Fluctuations and Pattern Formation in Self-Propelled Particles 
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We consider a coarse-grained description of a system of self-propelled particles given by hydrody- 
namic equations for the density and polarization fields. We find that the ordered moving or flocking 
state of the system is unstable to spatial fluctuations beyond a threshold set by the self-propulsion 
velocity of the individual units. In this region, the system organizes itself into an inhomogeneous 
state of well-defined propagating stripes of flocking particles interspersed with low density disor- 
dered regions. Further, we find that even in the regime where the homogeneous flocking state is 
stable, the system exhibits large fluctuations in both density and orientational order. We study the 
hydrodynamic equations analytically and numerically to characterize both regimes. 



I. INTRODUCTION 



Large collections of living organisms exhibit a highly 
coherent collective dynamics at large scales [lj . This be- 
havior, often referred to as "flocking" , spans an enormous 
range of length scales and is seen in diverse systems, in- 
cluding mammalian herds @ , crowds of pedestrians [1,3 , 
bird flocks ||, fish schools [f|, insect swarms Q, bac- 
terial suspensions Q, extracts of cytoskeletal filaments 
and molecular motor proteins 0, [uj , and motility assays 
pd| . Many of these systems can be unified under the 
theoretical paradigm of collections of self-propelled par- 
ticles. Their intriguing collective behavior has received 
considerable attention in recent years. 

A number of different theoretical approaches have 
proved fruitful in understanding the dynamics of collec- 
tions of self-propelled units. Starting with the seminal 
work of Vicsek [l2l |. rule-based models have been in- 
vestigated numerically and have been shown to exhibit 
noncquilibrium transitions between disordered and or- 
dered (flocking or moving) states. Subsequent work has 
focused on characterizing the nature of the order-disorder 
transition, its dependence on the noise, and pattern for- 
mation in the ordered state, both in the context of rule- 
based Vicsek-ty pe m odels [12h14| and of models of bac- 
terial swarming |27H2°| . Continuum hydrodynamic theo- 
ries have been used to describe the behavior of the system 
at large scales [l|, [HI, [IH, . Self-propelled particles are 
typically elongated and move along one direction of their 
long body axis. They can exhibit orientational order at 
high concentration. The ordered state is characterized by 
a vector order parameter, the polarization, which is also 
proportional to the mean velocity of the system. Hence 
the ordered state is a macroscopically moving state. The 
continuum theory has been developed phcnomenologi- 
cally on the basis of general symmetry arguments by 
drawing on analogies with magnetic systems and with 
liquid crystals [l|. In fact active or self-propelled sys- 
tems have been likened to "living liquid crystals" [la ]. 
This work has yielded several important results, includ- 
ing the possibility of long range order in 2D [16[ and 
the prediction and observation of giant number fluctua- 



tions in the ordered state [18|, |20j. The continuum theory 
has also been derived by systematic coarse-grainin g of 
specific microscopic models, including rule-based [2lll22l| 
and physically motivated [H, [3(| models. These deriva- 
tions yields (model-dependent) estimates for the param- 
eters in the hydrodynamic equations and have provided 
insight into the microscopic origin of of the large-scale 
collective physics. 

In a recent paper we derived the hydrodynamic equa- 
tions for a collections of self propelled hard rods mov- 
ing on a frictional substrate and interacting through ex- 
cluded volume interactions (23j . Although self-propulsion 
and steric effects alone are not sufficient to yield a ho- 
mogeneous polarized moving state in bulk, the hydro- 
dynamic equations are easily modified to incorporate a 
mean-field continuous transition from an isotropic state 
at low concentration of rods to a polar state at high den- 
sity. In the present paper we examine analytically and 
numerically the coupled nonlinear hydrodynamic equa- 
tions for density and polarization to characterize the 
large-scale structures that replace the linearly unstable 
homogeneous ordered state. The main results of our work 
are summarized in Fig. Q] that represents a "phase dia- 
gram" in terms of the density p of rods and their self- 
propulsion speed, vq. In the absence of self-propulsion 
(vq = 0) the model considered exhibits a mean-field 
continuous transition at the critical density p c from an 
isotropic state of zero polarization for po < p c to an 
ordered moving state, with uniform density and macro- 
scopic polarization P / 0. In a uniform ordered state 
at finite ^o, all rods would move with uniform mean ve- 
locity ~ voP. We find however that the moving state 
for po > Pc exhibits more complex behavior. For vq be- 
low a critical value v c (p ) the steady state of the sys- 
tem is still macroscopically polarized on average, but ex- 
hibits anomalous density and polarization fluctuations. 
We refer to this state as the "fluctuating flocking state" . 
The anomalous density fluctuations are the giant num- 
ber fluctuations predicted by Toner and coworkers [l7| 
and observed experimentally in active nematics |3lj . An 
additional feature of this regime, is a very slow tempo- 
ral approach to this noisy steady state, with some fea- 
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tures of a coarsening process. For vq > v c (po) the system 
orders in a robust striped phase, consisting of traveling 
high/low density stripes. The high density stripes are or- 
dered with polarization transverse to the long direction 
of the stripes, which is also the direction of the stripes' 
motion. Both of these phases have been identified in nu- 
merical studies of the Vicsek model [HI, [32|. Here we 
characterize them in terms of their origins in the hydro- 
dynamic equations using both analytical and numerical 
tools. This allows us to derive an understanding that 
transcends any specific microscopic model and is gener- 
ically applicable to a large number of self-propelled sys- 
tems. 

The coupling of density and polarization fluctuations 
embodied in the convective terms in the hydrodynamic 
equations of self-propelled systems plays a crucial role 
in controlling pattern formation. Some of these convec- 
tive terms reflect the dual role played by the polarization 
field as vector order parameter and as the mean velocity, 
resulting in competition between diffusion and convec- 
tion along the direction of mean local order. There is 
a qualitative analogy here with sedimentation problems 
[33l |34| , where the interplay of local alignment along the 
scdimcnting direction and diffusion can destabilize the 
system resulting in convective patterns, although hydro- 
dynamic interactions, not incorporated here, often also 
play an important role in sedimenting systems. 

The layout of the paper is as follows. First we intro- 
duce the hydrodynamic equations that are the starting 
point of our analysis. Then, we carry out a linear sta- 
bility analysis about the ordered state and characterize 
the region of linear stability of the bulk ordered phase 



The hydrodynamic equations (JTJ) and ([2j have the same 
form as those first proposed on a phenomenological basis 
by Toner and Tu [l|, Il6l . [T3| to describe the physics of 
flocking. The parameter ai is chosen to change sign at a 
characteristic density p c , while 04 > 0. This guarantees a 
mean-field continuous transition from an isotropic state 
with p = po and P=0 when 02 > to a homogeneous po- 
larized state with p = po and |Po| = \J when 02 < 0. 

These equations have also been derived from specific 
microscopic models of self-propelled particles on a sub- 
strate by some of us [23l - l25l | and by Bertin and collabora- 
tors [2l| . Bertin et al obtained hydrodynamic equations 
by coarse-graining a Vicsek-type model of self-propelled 
point particles, with a specific aligning rule for the pair 
interaction. In contrast, Baskaran and Marchetti, con- 
sidered a model of self-propelled hard rods of finite size 



or homogeneous flock. Next, we report the results of nu- 
merical solution of the nonlinear hydrodynamic equations 
and identify and characterize both the fluctuating flock- 
ing state and the striped phase, as well as the coarsening- 
like behavior leading to these phases. 

II. THE MODEL 

We consider a collection of polar rods of length £ mov- 
ing on an inert substrate characterized by a friction con- 
stant £ in two dimensions (2d) . Each rod is driven by an 
internal force F acting along one direction of its long axis, 
called its head. This force, together with the frictional 
interaction with the medium, results in a self-propulsion 
speed vq = F/( of constant magnitude. On length scales 
long compared to I and on time scales long compared 
to the microscopic interaction times, the dynamics of the 
system can be described in terms of hydrodynamic fields, 
namely the conserved densities (here the density p (r, t) 
of rods) and the variables associated with possible bro- 
ken symmetries. A collection of self-propelled polar rod 
can order in a polarized state, characterized by a finite 
value of a vector order parameter, P(r, t), describing the 
mean polarization of the rods. The ordered state is also 
a moving state, with mean velocity ~ vqP. The dynam- 
ics of the system is described by coupled equations for 
density and polarization, given by 

dtp = -V • (pv Q P - DVp) (1) 

and 



(2) 

I 

with excluded volume interaction and analyzed in detail 
the modifications induced by self-propulsion on the linear 
and angular momentum exchanged in a binary collision. 
The equations obtained by Bertin ct have precisely the 
form given in Eqs. ([T]) and ([2]), with parameters ai and 04 
determined by the aligning interaction between particles. 
In contrast, it was demonstrated in Ref. [23[ that steric 
effects alone are not sufficient to yield a homogeneous 
bulk polarized state. As a result, the equations derived 
in [23[ for a purely physical model have 04 = 0. We 
note that it was also demonstrated recently in a model 
of swimmers in a fluid that hydrodynamic interactions 
among swimmers are equally insufficient to yield a ho- 
mogeneous polarized state in bulk (26|. These results 
suggest that genetically and biochemically-regulated sig- 
naling, or external symmetry-breaking effects, such as 
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d t pP + \i P P ■ VpP = -D r [oa (p) + P 2 a 4 (p)] P P - ^-Vp + A 3 pPi VpP; + A 2 pPV • pP 
+ (D s - D b ) V (V • pP) + D b V 2 P P. 
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FIG. 1: (color online) Phase diagram in the (vo,po) plane. At vo — the system exhibits a continuous mean field transition 
at po — p c from an isotropic (I) to a homogeneous polarized (HP) state. The isotropic phase survives at finite vo in the region 
p < p c bounded by the vertical dashed line (blue online). For po > pc there is a critical v c (po) separating a polarized moving 
state with large anomalous fluctuations, named the "fluctuating flocking state" , at low self- propulsion speed from a high-speed 
phase of traveling stripes. The circles denote the values of v c (po) obtained numerically with the error bars indicating the step 
size used in the computation. The dashed-dotted line (purple online) is the longitudinal instability boundary Vci(po) obtained 
in Section III. The dashed line (red online) is the splay instability boundary v'c(po), given by Eq. (|A11|) . The top right panel 
is a real space snapshot of the density profile in the striped phase. The stripes travel in the direction of the white arrow, that 
also denotes the direction of mean polarization in the high density regions. The bottom right panel shows a real space snapshot 
of the density profile in the "coarsening" transient leading to the fluctuating flocking state at vo < v c . Density values from low 
(blue) to high (yellow) are indicated in the side bar. The value of the density in the blue stripes is well below the critical value 
p c — 0.5, while the red stripes are well into the polarized phase. 



chemotaxis, may be needed to obtain a polar state. Non- 
physical or external mechanisms of this type are embod- 
ied in Bertin et al in an alignment rule, but are absent 
in the work by Baskaran and Marchetti that aimed at 
identifying the role of purely physical interactions in con- 
trolling the large scale behavior of self-propelled systems. 
The goal of the present paper is to study the stability of 
the homogeneous polar state. For this reason we have 
added phcnomcnologically the term proportional to 
to the equations derived in 23). Both 02 and 04 will 
then treated as phenomenological parameters. 

The density satisfies a conservation law, with a flux 
controlled by two terms: pvqP describing convection 
along the mean self-propulsion velocity, vqP, and a dif- 



fusive current — Z?Vp that drives the system to a ho- 
mogeneous state. The anisotropy of the diffusion coeffi- 
cient relative to the direction of mean motion is neglected 
here for simplicity. The various terms in Eq. ^ (other 
than the one proportional to 04) are obtained from the 
microscopic hard rod model and have a simple physi- 
cal interpretation. The polarization field P plays a dual 
role in self-propelled systems. On one hand, it represents 
the vector order parameter associated with the sponta- 
neous breaking of rotational symmetry in the polarized 
state. Its dynamics is then in the class of that of equi- 
librium polar liquid crystals and X-Y spin systems. On 
the other hand, wqP is also the mean velocity of the flock 
with which particles are convected. The interplay be- 
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tween these two physical roles of the polarization field 
gives rise to the various terms in Eq. $2$ and underlies 
most of the phenomena discussed in this paper. The 
three terms proportional to A; in Eq. (J2|) play a crucial 
role in controlling the pattern formation phenomena de- 
scribed below. If we think of «oP as a velocity, then 
all three terms have the structure of convective nonlin- 
earities. Galilean invariance would require Ai = ^ and 
A2 = A 3 = 0. The self-propelled overdamped system 
considered here is, however, moving relative to a fixed 
substrate and docs not satisfy Galilean invariance. As 
a result, the values of Ai are unconstrained and in gen- 
eral model-dependent. There is an additional important 
difference between these three terms. The term propor- 
tional to Ai is a truly noncquilibrium term that can be 
understood only as a convective nonlincarity. In contrast, 
the terms proportional to A2 and A3 have an equilibrium- 
like interpretation associated with the role of P as the 
polar order parameter. These two terms would arise in 
equilibrium from a term of the form p |P| 2 V • P in the 
free energy, which effectively accounts for a dependence 
of the elastic constant associated with splay deformations 
on the amount of orientational order in the system. In 
this case one would obtain A2 = —A3. Finally, D s and D b 
arc diffusion constants that characterize the relaxation of 
splay and bend fluctuations, respectively, and D r is a ro- 
tational diffusion rate. 

Before proceeding to analyze the hydrodynamic equa- 
tions, it is useful to introduce dimensionlcss variables. 
We measure time in units of the inverse rotational diffu- 
sion rate, D^ 1 , and lengths in units of the length £ of the 
self propelled particles. The various dimensionless fields 
and parameters are then given by 

p = pl\ 
«b = v /(lD r ), 
A; = K/(D r £ 3 ), 



D = D/(D r £ 2 ). 

In the following all quantities are dimensionless and we 
drop the tilde for simplicity of notation. 



III. LINEAR STABILITY 

The hydrodynamic Eqs. ((T|) and ([2]) admit two homo- 
geneous solutions: an isotropic state (I) with p = po and 
P = for p < p c , and a homogeneous polarized (HP) 
state with p ~ po and P = -P0P0 for po > p c , where po 
is the direction of broken symmetry and Po = \J — a 2 j a\. 
The critical value p c is defined by 02 (p c ) = and is cho- 
sen here as p c — 0.5. It was shown in Ref. [23[ that the 
isotropic state is always linearly stable. In this section 
we examine the linear stability of the HP state at finite 
vq. To do this we linearize the hydrodynamic equations 
by letting 

p = p + 5p, (3) 
P = p (P + 5P) + P 5p±, (4) 

where po ■ Sp± = 0. Inserting this ansatz in Eqs. (QJ and 
([2|), we obtain three coupled equations for the fluctua- 
tions in the density, 6p, the magnitude SP of the polar 
order parameter and the director, Sp±. Combining the 
fluctuations into a vector, 

fSp(r,t)/p a \ 
Sy a (r,t)^\ SP(r,t) (5) 
V *Px(r,i) / 

and introducing the Fourier components, Sy (k, t) = 
J r e lk r (Sy (r, t), the coupled linear equations can be writ- 
ten in matrix form as 

5y a (k,t) = A a p(k)5y (k,t) , (6) 

where 



A(k) 



ik\\voPo — Dk 2 



-2aa 20 Pa + ik\\ (f + Apg P 2 ) 



-D<k 2 



-(D„-D b )knkj 



ikuVo 

2a 20 - ifc||ApoP 
- -D s fc|| 



-D b kl 



-ikx\ 3 p 



Pa 



ikj_v Q P 
-ik ± \ 2 p Po - (D a - P 6 )fc||fc_L 
ik\\Xip Q P Q - D s k\ - D b kf, 



(Z) 



with k = p fc|| + k^, k^ = k±k±, 5y 3 (k,t) = k± ■ where a 2 o = a2(po) < 0, a 40 = a 4 (/9 ) > 0, and 
<5pj_(k, t), and 

A = Ai-A 2 -A 3 . (9) 

a _ Po ( ^ a '2 \ _ Po ( 9&4 A ^ The coefficients a 2 and a 4 are chosen of the simplest form 

2<22o \ dp J p=po 2<i4o \ dp J p=po ' that guarantees a continuous transition at p c and Pq ~ 1 
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a 2 = 1 - p/pc 
a 4 = l + p/p c 



With this choice a 



POPa 

pI-p\ 



is always positive. 



for p > p c , i.e., 

(10a) 
(10b) 

,(k)t 

where the rates s a (k) are the hydrodynamic modes of 
the system. These are defined as those with decay rates 
(here proportional to Re[s a (k)]) that vanishes in the long 
wavelength limit k —> 0. Modes with Re[s a (k)] < 
decay at long times, while modes with Re[s a (k)] > 
grow, rendering the homogeneous state linearly unstable. 



We look for solutions of the form (5j/(k, t) ~ e t 



We discuss the hydrodynamic modes by considering some 
simplified cases. Further details are given in Appendix 
A. 

First, we consider the behavior of the system for 
Pa 3> Pc i-e., deep in the ordered state. The rate of 
decay of long wavelength fluctuations of the magnitude 
5P of the order parameter is controlled by A22 ~ 2a2o, 
which is always finite for p 3> p c , away from the mean 
field continuous transition. In other words SP is a non- 
hydrodynamic variable that decays on microscopic time 
scales. In this regime we can then neglect fluctuations SP 
and simply consider the dynamics of density and director 
fluctuations governed by the two coupled equations 



J 



d t 5p = (ik\\v Po - Dk 2 ) dp + ik±v p P 5pj 



(11) 



\2P 



PoPvk 



(D s - D b ) 



SJ 
Po 



r 



iknXxpoPo 



Spa 



(12) 



The general form of the dispersion relation of the hydro- 
dynamic modes is readily obtained by solving a quadratic 
equation and is given in Appendix A. Here we discuss 
some limiting cases. For wavevectors k along the direc- 
tion po of broken symmetry, i.e., k = k\\ and k± = 0, 
density and orientation fluctuations decouple and decay 
with rates 



s L p {k) =ikv P - Dk 2 
s% = ik\ lP oPo - D b k 2 



(13a) 
(13b) 



Both modes are stable and propagating, albeit with dif- 
ferent speeds. Deep in the ordered region, Pq ~ 1. The 



propagation speed of density fluctuations is then simply 
t>o, while the propagation speed of director fluctuations 
is Aipo ~ v oPo- A Galilean invariant system would have 
^i = vq/pq and the two modes would have the same 
propagation speed, vq. The difference in the propaga- 
tion speed of the two modes can then be considered a 
signature of the violation of Galilean invariance. 

Next, we consider wavevectors k transverse to the di- 
rection po of broken symmetry, i.e., k = k± and kit = 0. 
In this case the equations for density and director fluctu- 
ations are coupled and the two hydrodynamic modes are 
given by 



J 



"2^ 



D s )k 2 ± 



\{{D- D s f k 4 - 2k 2 v oPo [vo - 2 Pa P 2 X 3 ] } 



1/2 



(14) 



r 



The mode can become positive, yielding an instability, 
for k < k c , with 



y/vo [2p P^X 3 -v Q }/(DD s ) , (15) 



provided 



2P0-P0A3 > v 



The parameter A3 has been estimated for a 
croscopic models and found to be of order Vq 



(16) 
few mi- 

MM. 



Eq. (fT6j) then identifies a value (po) of the self propul- 
sion speed above which the homogeneous polarized state 



becomes unstable, with k c ~ (vq — uf ) 1 ^ 2 . The instabil- 
ity boundary i£ [po) depends on microscopic parameters 
and is therefore mo del- dependent. Using the parameter 
values obtained for the model of self-propelled hard rods 
discussed in [23j , where the nonlinear terms in the polar- 
ization equation arise from momentum-conserving colli- 
sions between the self-propelled rods, and summarized in 
Table HI we obtain 



[2irp P ( 



(17) 
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This instability is associated with splay deformations of 
the director and with spatial gradients normal to the di- 
rection of mean order, suggesting that it may bear no 
relevance to the stripe formation observed in the numer- 
ics due to the fact that the particles in the stripes are 
always aligned along the short direction thereby retain- 
ing the long-wavelength nature of the splay mode. How- 
ever, as shown in Fig. (fTJ), the splay instability line agrees 
remarkably well with the numerically observed onset of 
stripe formation at high density. This result, discussed 
further below, suggests that nonlinear pattern selection 
mechanisms may play an important role in stripe forma- 
tion. Finally, Eq. (fTTj) shows that ~ — for p >> 



the hydrodynamic mode takes the form 



Pc 



and appears to diverge as we approach the phase tran- 
sition. This apparent divergence is regularized when the 
effect of overdamped fluctuations of the magnitude of the 
polar order is incorporated in the mode analysis. This is 
done in Appendix A, where rather than just neglecting 
8P entirely, we approximate its behavior by assuming 
that on the time scales of interest |9t<5P| <C |2a2o^P| hi 
Eqs. ([6]) and (JT]). We then neglect the time derivative of 
5P, solve for the overdamped 8P. and use this result to 
eliminate it from the equations for density and director 
fluctuations. We find that the resulting modes still ex- 
hibit a splay instability as described above, but the crit- 
ical speed monotonically approaches the constant value 
v c(Pc) = [f /°c] _1 as Po Pc- The renormalized bound- 
ary uf(po) of this splay instability given in Eq. (|Alip is 
plotted in Fig. Q] as a red dashed line. In addition, fluctu- 
ations in the magnitude of the polarization renormalize 
of the diffusion constant associated with the decay of 
density fluctuations. 

As the continuous order-disorder transition is ap- 
proached from above, 020 — > and the separation of 
time scales between the decay of speed/magnitude fluc- 
tuations SP and the true hydrodynamic variables Sp and 
5p± no longer holds. To capture the physics of the system 
in the vicinity of the order disorder transition, we need to 
retain the dynamics of the "non- hydrodynamic" variable 
SP and examine the three coupled equations ((6J. One 
can show that the splay instability described above for 
k normal to the direction of mean order survives and is 
qualitatively unchanged. On the other hand, for k along 
the direction of broken symmetry, director fluctuations 
Sp± decouple from density and speed fluctuations and 
decay at the rate (|13bl) . The coupled modes for the dy- 
namics of density and speed fluctuations are then given 
by 



' ' n + A 22 )±-^/(A n - A 22 f + AA 21 A 12 (18) 



(A 



where Aij are the elements of the matrix A(k, t) given 
in Eq. for k = fciipo- It is easy to see that one of the 
two dispersion relations describes a non-hydrodynamic 
mode, s+ (0) = 2a 2 o, but with a decay rate that becomes 
vanishingly small for po — > p c . The other mode vanishes 
at k = 0. At small wave vectors the dispersion relation of 



s 1 : (k) = ikv Q P Q (a + 1) - D eff k 2 + 0(k 3 ) , (19) 



with 



D, 



eff 



D 



v 2 (a + 1) -kvIpq (a + 1) 



4|a 20 | 



2a 40 



2a 



in 



(20) 



where we have used the microscopic parameters given in 
Table 1. When D e ff < 0, density fluctuations grow in 
time and the ordered state is unstable. At the phase 
transition, i.e., for po = p c , The condition -D e // < is 
satisfied for all values of vq and the ordered state is always 
unstable. Away from the transition, by considering the 
exact modes in (|18p we find that for densities in a range 
Pc < Po < Pc there exists a range of self propulsion 
speeds v^i < «o < v c2 where the propagating density 
fluctuations are unstable. The lower instability boundary 
v ci(Po) i s shown in Fig. (TTJ) as a purple dashed-dottcd 
line. 

The results of the linear stability analysis are sum- 
marized in Fig. [5] The linear stability analysis predicts 
that near the mean field order-disorder transition the ho- 
mogeneous ordered state is destabilized at small vq by 
the growth of coupled density and polarization fluctua- 
tions. The instability occurs for spatial gradients along 
the direction of mean order, signaling the onset of spa- 
tial structures that are inhomogencous in this direction, 
like the stripes found numerically. The wavevector k c 
of the fastest growing mode for this instability scales as 
k c ~ (vq — Vc)^ 1 / 2 at fixed density po and as (po — Pc) 1 ^ 2 
at fixed self propulsion speed. The boundary of stabil- 
ity v c \(po) obtained from the linear theory vanishes at 
p c , in agreement with the onset of the striped phase 
obtained by numerical solution of the nonlinear equa- 
tions, as shown in Fig. (fTJ, but grows faster with vq 
than obtained numerically. This discrepancy is likely to 
stem from the fact that the full nonlinear dynamics of 
amplitude fluctuations must be incorporated to account 
for the behavior in these regions. In addition, the lin- 
ear stability analysis predicts that the homegeneous or- 
dered state is again stable at large self-propulsion speed 
for Vo > v^ 2 (po). This second line is shown in Fig. [2] 
The numerics, however, yield a striped phase in this re- 
gion. Finally, the longitudinal instability only exists for 
Po < Pc — 1-1 j while numerically stripes are observed 
at all densities above a critical velocity. Deep in the or- 
dered state, the linear stability analysis predicts that the 
homogeneous flocking state is destabilized by splay fluc- 
tuations of the order parameter. In this case the insta- 
bility is associated with spatial gradients in the direction 
normal to that the mean order. The wavevector of the 
fastest growing mode also scales as k c ~ (i'o — fc) 1 ^ 2 at 
fixed density po and as (po — Pc) 1 ^ 2 at fixed self propul- 
sion speed, but it is clear that nonlinear pattern selection 
mechanisms must be involved to yield the formation of 
the observed transverse stripes (associated with spatial 
gradients in the longitudinal direction) in this region. On 
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FIG. 2: (color online) The figure displays the linear stabil- 
ity boundaries in the (vo, po) plane. All lines have been cal- 
culated using the microscopic parameter values of Table [T] 
The vertical dotted line (blue online) is the mean field con- 
tinuous transition from the isotropic (I) to the homogeneous 
polarized (HP) state. The dashed-dotted lines (purple on- 
line) are the boundaries (calculated by numerical solution of 
D e f f(vo, po) = 0, with D e ff given by Eq. ([20jl) that define the 
region v^i < vo < v^2 where the homogeneous polarized state 
is unstable due to the growth of coupled density and polar- 
ization fluctuation associated with spatial gradients along the 
direction of mean order (longitudinal instability). The linear 
theory predicts that the homogeneous polar state is unstable 
in the ruled region to the right of the vertical mean-field tran- 
sition and bounded by these two lines. The dashed line (red 
online) is the splay instability boundary given in Eq. (IA11II . 
It terminates at a finite value at po — p c . The linear the- 
ory predicts that splay fluctuations destabilize the polar state 
in the entire ruled region above the dashed (red) line. The 
region where the system exhibits both the longitudinal and 
splay instabilities is cross-hatched. The longitudinal instabil- 
ity boundary v^\{po) vanishes for vo —t and po — > pt , in 
agreement with the numerics. 

the other hand, the instability line obtained from the lin- 
ear theory agrees remarkably well with the numerical on- 
set of stripes in this high density region (Fig. ([1])). More 
work is needed to understand stripe formation at high 
density and the origin of the associated length scale. 



IV. NONLINEAR REGIME 

To go beyond the linear stability analysis and investi- 
gate the nature of the flocking state above v c , we have 
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TABLE I: Diffusion constants and convective parameters for 
the model of self-propelled hard rods with excluded volume 
interactions discussed in Ref. [gjj. All diffusion constants are 
in units of £ 2 D r and all convective parameters are in units of 
£ 3 D r . The diffusion coefficients have been expressed in terms 
of the longitudinal diffusion constant Do of a long, thin rod. 
Below we use the low density value, Do = 1/4. 

solved numerically the full nonlinear hydrodynamic equa- 
tions. The numerical analysis has been carried out us- 
ing the specific parameter values obtained for the self- 
propelled hard rod model of Ref. (23[, summarized in 
Table |U All diffusion coefficients are enhanced by self 
propulsion of an additive contribution proportional to v 2 
that arises from the persistent nature of the random walk 
performed by Brownian, self-propelled rods. In the fol- 
lowing we discuss the properties of the system in terms of 
two dimensionlcss parameters, the self-propulsion speed, 
Vo, and the density of particles, po. In the numerics 
the coefficients d2 and 04 that control the continuous 
mean field phase transition from an isotropic to a polar 
state have been taken to be of the simple form given in 
Eqs. (jlObp with p c = 0.5 in units of the rod length. This 
form yields Po ~ (p — Pc) 1 ^ 2 for p p c and Po — > 1 for 

P » Pc- 

For generality we include fluctuations beyond the mean 
field level in the numerical analysis by adding Gaussian 
white noise terms in both the density and polarization 
equations of the forms V • f p (r,t) and fp(r, t), respec- 
tively. The random forces are chosen to have zero mean 
and correlations 

< f ip (r,t)f jp (r',t') >= %A p p(r,t)5(r - r')S(t - t'), 

(21) 

< f iP (r,t)f jP (v',t') >= ^-^^(r-r'Wi-O, (22) 

where A p and Ap are dimensionless noise strengths. The 
noise in the density equation scales as [^(r,^}] 1 / 2 , while 
the polarization noise scales as [p{r,t)]-V 2 g||3||. This 
difference arises because the fields p and P are extensive 
and intensive quantities, respectively. The numerical re- 
sults described below are all for fixed values of the noise 
amplitudes, A p = A p — 0.3. We have solved the non- 
linear equations using the Euler method for numerical 
differentiation on a grid with Ax = 1.0 and At = 0.1 (we 
have verified that the numerical scheme is convergent and 
stable for At /(Ax) 2 < 0.5) We consider a square system 
of size L x L with both periodic and shifted boundary 
conditions and a range of system sizes. 

The behavior of the system as a function of the self- 
propulsion velocity vq and the density of particles p is 
summarized in the phase diagram shown in Fig [1] dis- 
cussed in the Introduction. The isotropic state is stable 
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for all t>o and p < p c . For p > p c and vq < f c (po) 
the system is in the fluctuating flocking state, charac- 
terized by finite polarization and large spatial and tem- 
poral fluctuations of both density and order parameter. 
For v > v c (po) we find a striped phase, with alternating 
ordered high density bands and disordered low density 
bands, propagating in the direction of order. In the nu- 
merics the value of v c is identified as the self propulsion 
velocity where the density histograms shown in Fig. (|9]) 
change from unimodal to bimodal. The histograms are 
constructed by recording the local density at each spatial 
grid point for fixed mean density po averaged over many 
initial conditions. We have also verified that histogram 
of local polarization magnitude change from unimodal to 
bimodal as the same value of vo- The numerical bound- 
ary for the onset of the stripe regime vanishes with vq for 
Po Pc and is in qualitative agreement with the bound- 
ary calculated in section for the onset of the longitudinal 
instability. The theoretical curve, however, grows much 
faster with vq than the numerical boundary. Surprisingly, 
at high density the theoretical boundary for the linear in- 
stability of splay fluctuations agrees very well with the 
numerical onset of stripes. 




A. Fluctuating nocking state 

In this subsection we characterize the fluctuating flock- 
ing state that exists in the region p > p c and Vq < v c 
of the phase diagram in Fig. [1] As noted earlier, this 
state is characterized by large fluctuations in the den- 
sity. These fluctuations do not, however, destroy the 
underlying orientational order of the system. This is dis- 
played in Figs. [3] and @] Figure |3] shows the time evo- 
lution of the magnitude squared of the order parameter, 
(P 2 (t)) = (P%(t) + Py(t)}, where here and below the 
brackets denote a spatial average over the system and an 
average over different realizations of initial conditions, 
for both an initial ordered {{P 2 {t = 0)) = 1) and an 
initial disordered ((P 2 (t = 0)) = 0) state. Both states 
reach the same ordered state at long times, although on 
very different time scales. The asymptotic state is or- 
dered (i.e., (P 2 ) ^ 0)and this does not appear to be 
an artifact of the finite system size, as shown in Fig. |4] 
where the magnitude of the order parameter is displayed 
as function of system size for various values of i>o and dif- 
ferent initial states. The numerics suggest that P is finite 
and close to unity (note the narrow range of (P) on the 
vertical axis of Fig. |4} in the fluctuating flocking state. 
The large difference in the relaxation time from an ini- 
tial disordered/ordered state to the asymptotic ordered 
state seen in Fig. [3] is not unexpected. When starting 
in a disordered state, the system locally finds different 
degenerate ordered states, which subsequently coarsen 
towards the homogeneous ordered state. Below we char- 
acterize both the coarsening behavior as the system seeks 
out its asymptotic steady state and the properties of the 
asymptotic fluctuating flocking state. 



FIG. 3: (color online) The magnitude squared of the orien- 
tational order parameter (P 2 (t)} as a function of time for an 
isotropic initial state ({P 2 (t — 0)} =0) and an ordered initial 
state ((P 2 (t — 0)) = 1), for three system sizes and vo = 0.1, 
po = 0.7, corresponding to v c ~ 0.42. The three curves ob- 
tained for an isotropic initial state overlap and cannot be 
distinguished in the figure. Both initial states approach the 
same macroscopically ordered state at long time, although the 
time scale required for the isotropic initial state to reach the 
asymptotic steady value is much longer. 

To quantify the coarsening behavior we have measured 
the two-point correlation function of both the density and 
the order parameter, defined as 

C p (r, t) =< Sp(v + r, t)5p(r , t) > , (23) 
C P (r, t) =< P(r + r, t) • P(r , t) > . (24) 

Before discussing the behavior of these correlation func- 
tions, it is useful to recall the dynamics of phase order- 
ing developed in the context of equilibrium second order 
phase transitions (37j . Phase ordering theories consider 
a system in an initially disordered state that is rapidly 
quenched below the order-disorder transition point and 
describe the time evolution following the quench. Imme- 
diately after the quench the system consists of finite-size 
ordered regions, each in one of the continuum of degen- 
erate ground states that correspond to one choice of the 
spontaneously broken continuous symmetry. The sys- 
tem then evolves in time and "coarsens" with some of 
the ordered regions growing at the expense of others and 
eventually taking over the entire system. The coarsening 
process is typically controlled by a single energy scale, 
namely the energy cost of the domain walls between dif- 
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FIG. 4: (color online) The mean magnitude of the polar order 
parameter as a function of system size L for different values 
of vq and p = 0.7, corresponding to v c ~ 0.42. The curves 
appear to saturate as the system size increases, suggesting 
that the system remains macroscopically ordered in the ther- 
modynamic limit. 



ferent ordered regions. This implies that the time evo- 
lution of the system occurs via the growth of a single 
length scale L{t) that characterizes the size of a typical 
ordered region in the system at a time t. The order pa- 
rameter correlation function will then depend on time 

only through L(t), i.e., Cp(r,t) = C^-^y^. Scaling 

analysis indicates that the length scale L(t) grows with 
the dynamical critical exponent z, L (t) ~ t x l z . For a 
vector order parameter in two dimensions in equilibrium 
one expects z ~ 2 [371 ], 

In the case of self propelled particles, the orientational 
fluctuations that drive the coarsening of the ordered state 
also induce mass fluxes and hence couple to density fluc- 
tuations. As a result, density correlations are essentially 
slaved to the order parameter correlations and both C p 
and Cp are expected to exhibit coarsening behavior [38j]. 
This is indeed the behavior that has been observed in 
active nematic liquid crystals, where both density and 
orientational correlations have been shown to coarsen on 
a characteristic length scale that grows like t 1 ^, with 

2~ 2 d|. 

The behavior of polar active system appears to be 
somewhat different. Figure [5] shows the two-point cor- 
relation function for the density and the order parameter 
for our flock. Both exhibit a growing correlation length 
as a function of time, but they do not exhibit the simple 



scaling behavior outlined above, indicating that the ap- 
proach to the homogeneous state is no longer controlled 
by the single energy scale associated with the cost of a do- 
main wall. Convective fluxes induced by self-propulsion 
lead to correlations on longer length scales and hence ac- 
celerate the coarsening dynamics. This picture can be 
substantiated by extracting a length scale L (t) from the 
correlation functions. This is shown in Fig. [S] where we 
see that the dynamical exponent z is smaller than the 
equilibrium value, indicating that the coarsening dynam- 
ics in polar active systems is faster than that of both 
equilibrium systems with vector order parameters and 
active nematics. This may be due to the fact that, in a 
polar self- prop elled system, there is true long range or- 
der in 2D [l6l. [l7j and the associated suppression of the 
Goldstone mode by the nonlinear couplings changes the 
dynamics of the system as the orientational order builds 
up towards the homogeneous ordered state. 

At long times, for vq < v Cl the system reaches the 
"fluctuating flocking state", a steady state with finite 
mean polarization and anomalous fluctuations. To char- 
acterize the properties of this state we have evaluated 
the two-point correlation function of fluctuations in the 
orientational order parameter, 



C SP (r, t) = < <SP(r + r, t) • <5P(r , t) >, 



(25) 



with <5P(r, t) = P(r, t) - P(t). This is shown in the right 
frame of Fig. [8] and it decays logarithmically as expected 
for vector order in 2D. The correlation function shown 
in Fig. [5] has been averaged over rrj, hence it represents 
only the isotropic (angular-averaged) part of the order 
parameter correlations. In general we expect the corre- 
lation function to be anisotropic and its spatial decay to 
be described by different length scales in the directions 
longitudinal and transverse to the direction of mean mo- 
tion, as described in Ref. [l7j. We have calculated the 
spatial decay of the order parameter correlation in each 
direction and find it indeed to be anisotropic, as shown 
in Fig. (JT)). The theoretical analysis of Ref. [131 predicts 
a power law behavior with different exponents character- 
izing the decay along and orthogonal to the direction of 
broken symmetry. Given the large spatial and temporal 
fluctuations in our system, the system sizes considered 
here are too small to obtain reliable statistics to quantify 
this behavior and extract scaling exponents. 

In contrast to the order parameter correlations, which 
decay logarithmically, the two point density correlation 
shown in the in the left frame of Fig. [5] displays correla- 
tions over longer length scales than expected in equilib- 
rium. In fact the density correlation functions exhibits 

cuspy behavior of the form C p (r, t) = 1 — ^ jpyl with 

a ~ 0.6 typically characteristic of a state with growing 
domains. Furthermore, C p (r, t) depends very weakly on 
the self propulsion speed. These large correlations in the 
density arise because of the coupling of this conserved 
field to the fluctuations in the underlying order parame- 
ter field that is an intensive variable. This leads to what 
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FIG. 5: (color online) Early time, two-point density (top panel) and order parameter (bottom panel) correlation function for 
vo = 0.1 < v c — 0.42 mean density po = 0.7, and L = 1024. The dashed horizontal line indicates the value of the correlation 
function at which we extract the coarsening length scale L(t), shown in Fig. [6] 
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FIG. 6: (color online) The coarsening length L(t) as a function 
of time for the density (circles) and order parameter (squares) 
correlations for Vo = 0.1, po = 0.7 and L — 1024. For this 
density v c — 0.42 and the system was started in a disordered 
state. The straight line has slope 2. Although no single scal- 
ing exponents can be extracted for the two length scales L(t), 
it is clear that the growth with time is faster than would be 
obtained for z — 2. 



has been termed "giant number fluctuations" in these ac- 
tive systems [20l | and is the underlying mechanism for the 
large fluctuations in the density in the flocking state of 
our system. 



B. Striped Phase 

The top inset of Fig. Q] shows a real space snapshot of 
the density obtained for po > Pc and vo > v c (po) In this 
region the systems consists of well defined stripes of the 
high density ordered phase alternating with stripes of the 
low density disordered phase. In the ordered region the 
polarization is always normal to the long direction of the 
stripes and the stripes travel at a fixed speed in the di- 
rection of polar order. Panels (a) and (b) of Fig. [9] show 
histograms of the density and magnitude of the order pa- 
rameter for a fixed density pa = 0.7 > p c and different 
self propulsion speeds, vq. For small Vq, the histograms 
are unimodal (fitted with a gaussian peaked at mean den- 
sity po = 0.7), signalling a uniform state. Above a char- 
acteristic value of vq the histograms acquire a bimodal 
structure (fitted with two overlapping Gaussians peaked 
at low and high densities), corresponding to the striped 
phase. The boundary v c (po) corresponding to the onset 
of the striped phase and shown in Fig. [1] is determined as 
the value of vo corresponding to the onset of this bimodal 
structure. The error bars on these data point are sim- 
ply the step size of our increments in vq. Within these 
error bars, the same values of V c (po) are obtained from 
the onset of a bimodal structure in both the density and 
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FIG. 7: (color online) Plot of Csp (r, t) for the case when r is 
parallel to po (black/open circles) and when r is perpendicular 
to it (red/ filled circles). The anisotropy of the correlation 
functions begins to appear at large r, but larger system sizes 
are needed to quantify the difference. 

polarization histograms. The points shown in the phase 
diagram are obtained using the density histograms. As 
noted earlier, this boundary closely tracks the thresh- 
old line for the onset of the linear instability of splay 
fluctuations at large po (although this instability arises 
from spatial gradients normal to those of stripe forma- 
tion!), and vanishes with vq as po — > pt > as predicted 
by the longitudinal linear instability discussed in Section 
III. However, there is considerable discrepancy between 
the linear theory and the numerics in the details of the 
behavior at low density. 

The bottom two panels of Fig. |9] show histograms of 
density and polarization for a fixed self propulsion speed 
v > v c and three different values of density. These his- 
tograms are used to infer the properties of this striped 
phase. The difference in position of the two peaks in 
the bimodal histograms indicates the contrast in density 
and order parameter between the ordered and disordered 
stripes. The position of the peaks is independent of the 
self propulsion speed and only weakly dependent on the 
density, suggesting that the density contrast between the 
isotropic and ordered stripes is entirely diffusion limited. 
The height of the high density/high order peaks in the bi- 
modal histograms is a measure of the width of the ordered 
stripe with respect to the disordered one. We note from 
the figure that the height decreases with both increasing 
vq and po- This indicates that the width of the ordered 
stripe decreases with increasing values of these two pa- 
rameters. Further, we measure the speed of propagation 



of the stripes, shown in Fig.QJJ As naively expected, the 
propagation speed of the stripes increases linearly with 
v . 

An alternative way of displaying the existence of the 
stripes and quantifying their properties is provided by the 
two point density correlation function defined in Eq. (|23|l . 
We have evolved the system starting form a uniform or- 
dered initial state and evaluated the two-point correla- 
tion function as a function of r for three directions: 0°, 
45° and 90° to the direction of initial orientational order. 
The result is shown in Fig. [TT] When the self propul- 
sion speed vq < v c (right frame), the correlation decays 
monotonically in all directions. On the other hand, when 
> v c (left frame), we find well defined oscillations in 
the correlation function in the directions normal and at 
45° to the direction of motion showing the emergence of 
the periodic structure associated with the stripe pattern. 
Also, this clearly indicates that the spatial inhomogencity 
develops in the direction of initial orientational ordering 
even in the region where the linear instability is along a 
wavevector orthogonal to the ordering direction. Finally 
we have also investigated the dependence of stripe for- 
mation on system size and boundary conditions. We find 
that the width and speed of the stripes remain mostly 
invariant as we go to larger system sizes. We have also 
solved our equations with shifted boundary conditions 
[4fj| and have found that the striped phase persists. 

In summary, for values of vq > v c (p) given by the 
(black) solid line in Fig. (JJJ), the system develops robust 
propagating stripes of alternating ordered and disordered 
regions. The numerically identified transition line follows 
closely the threshold for the splay instability identified in 
Eq. ([TT)) . goes to zero at the phase transition in agreement 
with the longitudinal instability identified using Eq. ([20)1 
and shows a behavior unlike both of these linear insta- 
bilities in the intermediate region. Further, for systems 
initialized in a uniform ordered state, these stripes form 
along the direction of initial ordering even in the domain 
where the longitudinal instability is absent. This sug- 
gests that the pattern selection arises from a complex 
interplay between the unstable linear modes [4l[. Also, 
the width of the stripes exhibits a scaling behavior consis- 
tent with the critical wavevectors k c of both instabilities, 
but is not quantitatively captured by either length scale. 
A systematic study of the relationship between the linear 
modes and the patterns observed here will be the focus 
of a future work. 

We can show that the nonlinear equations admit a 
propagating front solution that may correspond to the 
onset of the stripe phase, although the stability if this 
solution is yet to be established. In our numerical study 
we have solved the equations (UJ and © by systemati- 
cally dropping various nonlinear terms. We have estab- 
lished that the terms that are critical for the formation 
of the striped phase are the homogeneous nonlinearity 
in the coefficient that induces the phase transition, 
the couplings between density and polarization embodied 
by the convective terms in the density and polarization 
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FIG. 8: (color online) Late-time two point correlation functions of density, C p (r), and order parameter, Csp(r), fluctuations. 
The data are obtained from an initial ordered configuration for po = 0.7, L = 512 and various values of vo, in the region 
corresponding to the fluctuating flocking state. Inset (top panel) shows (1 — C p (r/ L{t))) vs. the scaled distance r/L(t), for 
various values of vq. The dashed straight line shows the cuspy nature of two-point correlation function with cusp exponent 
a = 0.6. Inset (bottom panel) shows the scaled two-point order parameter correlation function, Csp(r/L(t)), vs. the scaled 
distance r/L(t), for several vq. The dashed straight line shows the logarithmic decay of Csp(r / L(t)). 



equations (i>oVpP and v (Vp)/p, respectively) and the 
convective nonlinearity controlled by the parameter A3. 
The longitudinal instability arises from the interplay of 
(Z2 and the convective terms while the splay instability 
is controlled by the term proportional to A3. It is useful 
then to consider a simplified description of the nonlin- 
ear dynamics where diffusion is neglected and only terms 
essential for the pattern formation are retained, given by 



d t p = - V ■ (pv P) 
d t P = - [02 (p) + P 2 a 4 (p)] P 



2p 



Vp + XsPiVpPi 



(26) 



(27) 



Denoting by x the direction of broken symmetry of the 
putative HP state, we postulate a solution of these equa- 
tions in the form of a front uniform in y and propagating 
along x with a yet undetermined constant speed U, 

p(x,y,t)=p(x-Ut), (28) 
P(x,y,t)=P(x-Ut)x. (29) 

Inserting this ansatz, Eqs. (|26|) and (|27|) become 



d x \np 



(U - v P) 



d x P, 



(30) 



''0 



2p 



a 2 + ai P z ) P- A3P - d xP -(U + X 3P P) d x P = 0. 



(31) 



The density equation can be formally integrated by pos- 
tulating an isotropic state at x = 00 (P(oo,t) = 0) and 
a polar state at x = —00 (P(— 00, t) = 1). This gives us 
the ratio of the density in in a polarized region to the 
density an isotropic region as 



Ppol _ U 
Pi so U - V P 



(32) 



The density contrast is infinitely sharp when the front 
propagates at a speed vqP commensurate with the de- 
gree of ordering in the stripes. Diffusion, neglected here, 
will smooth the density crossover between the two re- 
gions. This is in agreement with our numerical results 
that indicated that the contrast between the two regions 
is insensitive to the parameters and is indeed diffusion 
limited. 

Writing the density from Eq. (|30|) as p = 1/(U — vqP), 
and substituting in the order parameter equation, we can 
formally integrate Eq. (f3~T|) to obtain a solution of the 
form 



A 



Aa4 



(P-l) 



In [ a4 (P-l) 2 (l + P)]+-^ln(P), 
Cr 2a 2 

(33) 

where for simplicity we have assumed U ~ vq and we 
have introduced a dimensionless friction constant, Cr = 
(02 + ^4) and a dimensionless length scale, A = . 
This formal solution cannot be inverted analytically. A 
plot is shown in Fig. Q2] and clearly displays that the 
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FIG. 9: (color online) (a)-(c) Density and (b)-(d) order parameter histograms. In panels (a) and (b) the histograms are shown 
for pa = 0.7, Vq = 0.1, 0.3, 0.5, 0.7 and L = 128. For small vo, below the threshold value v c ~ 0.42, the histograms are unimodal 
(fitted by a Gaussian), indicating a uniform state. For vo > v c , the histograms are bimodal (fitted by two Gaussian curves), 
indicating the onset of stripes. In panels (c) and (d) the histograms are shown for vo > v c and three values of the mean density, 
po = 0.6, 0.7, 0.8. The position of the peaks does not change with density, while the difference in the height of the two peaks 
increases with increasing density. 



solution represents a propagating domain boundary of 
effective thickness A between a state with P = 1 and 
a state with P = 0. In physical units, the length scale 
controlling the crossover between isotropic and polarized 
states, hence the sharpness of the stripes, is given by 
(3wvo/£D r ). This length is essentially the distance trav- 
eled by a self propelled particle in a rotational diffusion 
time. In other words, stripe formation is controlled by 
the formation of domain boundaries in the polarization, 
and the fact that the density is slaved to the polarization 
and hence leads to mass fluxes that delineate the two re- 
gions in the striped phase. Finally, it is apparent from 



Eqs. ([30|l and (j3l"j) that there is no propagating front so- 
lution if we turn off the couplings between density and 
polarization. 



V. SUMMARY 

In this work we have considered a continuum descrip- 
tion of a collection of self propelled particles moving in 
a passive medium. Their dynamics on large length and 
time scales is governed by hydrodynamic equations for 
the density and the polarization field. The crucial physics 
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FIG. 10: (color online) Speed of the stripes as a function of 
Vo — v c for po = 0.7. The speed of propagation of the stripes 
increases linearly with vq- The solid line is a guide to the eye. 
The dashed line is a linear fit ~ vq — v c 



in these systems that distinguishes them from conven- 
tional liquid crystalline systems is the dual role of the 
polarization field as i) a physical velocity that leads to 
mass convection and hence couples orientational fluctu- 
ation to density fluctuations and ii) an order parameter 
associated with a spontaneously broken continuous sym- 
metry. This duality leads to the remarkable phenomenon 
of long range ordering identified in [16]. Here we show 
that this same physics destabilizes the homogeneous or- 
dered state above a critical value of self propulsion speed 
and allows the nonlinear equations to admit a propagat- 
ing front solution that yields the striped phase identified 
numerically. Further, the coupling of orientational fluc- 
tuations to density fluctuations gives rise to anomalous 
fluctuations even in the regime where the ordered state is 
stable and leads to nontrivial coarsening dynamics that 
is different from the dynamics of both the equilibrium 2D 
X-Y model and that of active nematics. 

The two phases observed here, namely the striped 
phase and the fluctuating flocking phase, have been iden- 
tified earlier in the context of numerical studies of the 
Vicsek model. Our work identifies the origin of these 
phenomena in the model independent framework of the 
dynamics of conserved quantities and broken symmetry 
variable. It has been shown in different systems of this 
class that pattern formation phenomena might be cru- 
cially related to biological functionality [42|, |43|. This 
work would facilitate the application of theoretical tools, 
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FIG. 11: (color online) Steady-state two-point density cor- 
relation function for three different directions with respect 
to the ordering direction, 90°, 45° and 0°, for va > v c (top 
panel) and for vq < v c (bottom panel). For vq < v c , there 
is no directional dependance in the correlation function. For 
Vo > v c , correlations at 90° to the ordering direction decay 
monotonically, while correlations at 0° and 45° to the order- 
ing direction show oscillations. 



such as the amplitude equations and pattern selection 
analysis that are well developed in the context of chem- 
ical reacting systems to collections of self propelled par- 
ticles. 



15 



discussions. 




FIG. 12: Plot of polarization as a function of x as obtained 
by inverting Eq. [33] for po = 0.7 and vo = 1.0. The solution 
represents a propagating domain boundary between a state 
with P — 1 and a state with P = 0. For these parameters 
A ~ 1.57. 



Appendix A: Linearized Equations 

Here we analyze in more detail the hydrodynamic 
equations linearized about the homogeneous polar state 
given in Eq. ([6]) and (0 in the main text and present 
a better approximation for the discussion of the splay 
instability. 

The hydrodynamic modes deep in the ordered phase 
were discussed in Section III by entirely neglecting mag- 
nitude fluctuations SP that decay on microscopic time 
scales. A better approximation consists of neglecting the 
rate of change of SP in Eq. (??) and solve for SP in 
terms of fluctuations in density and director to lowest 
order in gradients, with the result 



SP = 



Sp 
Po 



(Al) 



We then use this expression to eliminate SP from from 
Eqs. (|6]) for density and director fluctuations. The eigen- 
values of the resulting two coupled equations are given 
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^(cii + c 22) ± 2\A c n ~ C 22) 



4ci 2 C2i (A2) 



J 



di = ik\\v P (1 + a) - fen 



D 



1 v, 



2a 2 o 



(A3) 



C12 = ik^VoPoPo - fc||fc_LPo w O- P TT" 2 - 

Zd 



20 



(A4) 



c 2 i = ik± 



- P0A3P0 (1 + a) 

ZP Q 



- kuk± 



(D s - D b ) - 



A3P0 (Vq 



2a 



■2(1 



(A5) 



C22 = -ifc||Ai/o Po 



D b k\ 



n _ ^2Pa p2 \ ,2 
2a2o 



(A6) 



Again the modes decouple when k = fciipo lies along the the dynamics of density and director fluctuations are then 
direction of broken symmetry. The two modes governing 
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given by 



ikv P (1 + a) 



D 



4a 



2ii 



2Xv p P 



= -ik\xp Pu - D b k 



k 2 , 
(A7a) 
(A7b) 



Director fluctuations are stable and their decay rate is 
controlled by the bend diffusion constant. Since 0,20 < 
in the ordered state, one can define an effective dif- 
fusion constant in Eq. (|A7a[) as D^ff = D + (Vq + 
XvoPqPq ) / (4|a2o I ) ■ The first correction to D in this ex- 
pression, proportional to v 2 , always enhances diffusion 
and arises from the fact that self propelled particles per- 
form a persistent random walk |23(. The second correc- 
tion can lead to an instability if A < 0. The parameters 
Ai are microscopic quantities and their values are model 
dependent. As discussed in the main text, if we think 
of the polarization an an equilibrium order parameter, 



then A = 0. If in contrast we think of the polarization as 
a physical velocity in a Galilean invariant system, then 
A > 0. In both these cases the density fluctuations relax 
diffusively for all values of the parameters. For the self- 
propelled hard rod model discussed in Ref. [23| A > (see 
also Table HJ. In this case the convective terms propor- 
tional to Ai further enhance the effective diffusion con- 
stant and the homogeneous state is stable. Note that a 
value A > is also obtained in the microscopic Boltz- 
mann equation model studied in pH |. If, however, the 
microscopic model allows for higher order chemical and 
biological processes that can lead, for example, to a re- 
versal of the direction of motion of an individual unit due 
to interactions with other units, then A can be negative 
and drive a longitudinal instability [35j . 

Next, we consider wavevectors k transverse to the di- 
rection po of broken symmetry, i.e., k = k±. In this case 
the equations for density and director fluctuations are 
coupled and the two hydrodynamic modes are given by 



J 



"D + D.) k 2 ± i{ (D - D s ) 2 k 4 - 2k 2 v oPo [v - 2 Po P 2 X 3 (1 + a)] } 



1/2 



(A8) 



where D s = D s + ^P 2 A 2 A 3 /(2|a2o|) 
become positive, yielding an instability, for k < k, 



The mode sT can 
, with 



provided 



v [2p p2\ 3 (l + a)-v }/(DD s ), (A9) 



2p PoA 3 [l + a] > v . (A10) 



Using the parameter values obtained for the model of self- 
propelled hard rods discussed in p3j , where the nonlinear 
terms in the polarization equation arise from momentum- 
conserving collisions between the self-propelled rods, and 
summarized in Table HI we obtain 



vf = [2Trp P 2 (l + a)Y 



(AH) 



v c (Pc)- Finally, for a wavevector k at an angle 9 to direc- 
tion po of broken symmetry, the splay instability occurs 
for a range of angles 9 m < 9 < tt/2, where 9 = ir/2 
corresponds to k normal to po- The growth rate of the 
unstable mode is, however, always largest for 9 = tt/2, 
when director deformations are pure splay. 
Appendix B: Wavevector of fastest growing modes 



In this appendix we identify and charachterize the 
fastest growing mode associated with the two linear in- 
stabilities identified in the main text. To identify the 
wavevector of the fastest growing mode for the longitu- 
dinal instability discussed in Section III, we evaluate the 
real part of the dispersion relation of this mode to order 
fc 4 , with the result 



The critical line (po) given in Eq. (|ATTj) is plotted 111 
Fig. HI As obtained in the main text , the instability line 
vanishes as vf ~ 1/po at large density. However, near the 
transition incorporating overdamped magnitude fluctua- 
tions regularizes the behavior yielding a finite value for 



Rels^k)} = -~D eff k 2 - D 4 k 4 + 0(k 6 ) , (Bl) 



where D e f t is given in Eq. (f2T)f and 
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D 4 = 



1 



32 | a 2 o 
1 



^A{A-4aP v ){D-D s ) 
1 



32 |a2o|' 
3 



8 |a 20 



16 8|a 20 | 



{A-4aP v Q y 



A 1 -2 (D - D s ) - 



|«20 
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8 |a 20 
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FIG. 13: (color online) The maximum growth wavevector of 
the longitudinal instability as a function of the self-propulsion 
speed Vo for different densities. 



opposite behavior. This indicates that the length scale 
selected by the nonlinear pattern is not simply related 
to the wavevector of the most unstable mode associated 
with the linear instability. 

kc 

25 



0.20 



0.15 



0.05- | 




where 

A = 2[v (l-2a)+\p ]P (B3) 

The growth rate of the unstable mode is then maximum 
at a wavevector k c , given by 

k L c = yJ-D eff /2D 4 (B4) 

Fig. (|13[) shows a plot of the maximum growth wavevec- 
tor as a function of the self propulsion speed for vari- 
ous values of the mean density po. The critical length 
scale k~ x decreases with increasing density poj m agree- 
ment with what observed numerically for the width of the 
stripes. On the other hand, fc" 1 increases with increas- 
ing ^decreases as the self propulsion speed increases, im- 
plying that the width of the stripe should increase with 
increasing SP speed, while the stripes width exhibits the 



FIG. 14: (color online) The maximum growth wavevector of 
the splay instability as a function of the self-propulsion speed 
Vo for different densities. 



Next, proceeding as above, we can find the fastest 
growing mode associated with the splay instability. This 
is of the form 

(D s + D) V \v c J { ' 

This critical wavevector is shown function of self- 
propulsion speed vo > t>f (p) for different values of po in 
Fig. (|14[) . In this case k c is a non-monotonic function of 
Vq. But, it increases with both Vq and po for the range 
of densities and self propulsion speeds probed by the nu- 
merical analysis. 
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